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Abstract 

We present some applications of an Interacting Particle System 
(IPS) methodology to the field of Molecular Dynamics. This IPS 
method allows several simulations of a switched random process to 
keep closer to equilibrium at each time, thanks to a selection mech- 
anism based on the relative virtual work induced on the system. It 
is therefore an efficient improvement of usual non-equilibrium simula- 
tions, which can be used to compute canonical averages, free energy 
differences, and typical transitions paths. 

Keywords: Non-equilibrium molecular dynamics. Interacting Particle 
System, Genetic Algorithms, Free energy Estimation. 
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Phase-space integrals are widely used in Statistical Physics to relate the 
macroscopic properties of a system to the elementary phenomenon at the mi- 
croscopic scale il4j . In constant temperature (NVT) molecular simulations, 



these integrals often take the form 

^^{A) = {A)= [ A{q,p)dfi{q,p). (1) 
Jt*m 

where M denotes the position space (also called the configuration space), 
and T*M denotes its cotangent space. A generic element of the position 
space M. will be denoted by g = {qi,- " , Qn) and a generic element of the 
momentum space hy p = {pi, ■ ■ ■ ,pn)- We will consider here that M ~ 
or (a torus of dimension 3A^, which arises when using periodic boundary 
conditions), and that T*M ~ M^^ x R^^ or T^^ x M^^, though in general 
more complicated situations should be considered, when performing Blue 
Moon sampling [SI, IS] for example. 

The measure // is the canonical probability measure 

dfJ-{q, p) = ex.p{-(3H{q, p)) dq dp, (2) 

where (3 = l/k^T (T denotes the temperature and the Boltzmann con- 
stant) and where H denotes the Hamiltonian of the molecular system: 

Hiq,p) = ^p''M-'p + V{q). (3) 

In the above expression, V is the potential experienced by the particles, 
and M = Diag(mi, • • • , m^v) where rrii is the mass of the i-th particle. The 
constant Z in @ is the normalization constant defined as 

Z= ex.p{—PH{q,p)) dqdp. 

Jt*m 

Some quantities can not be expressed through relations such as One 
important example is the free energy of a system, defined as 

F = -13-^ In Z. (4) 

It is though often the case in practice that a straightforward sampling of 
IJL is difficult. Indeed, high dimensional systems exhibit many local minima in 
which the system remains trapped, especially when the temperature is low. 
In those cases, alternative approaches have to be used, such as those built 
on the simulated annealing [23[ paradigm. The idea is to switch slowly from 
an initial simple sampling problem, to the target sampling problem, through 
a well chosen interpolation. This allows to attain deeper local minima, but, 
due to its nonequilibrium nature, is not efficient as such to sample accurately 
the target measure. 

We mention that variations have been proposed, especially tempering 
methods (see |19| for a review), the most famous being parallel temper- 
ing [2Z|, which also has an importance sampling version [2H]- These meth- 
ods consider an additional parameter describing the configuration system 



2 



(e.g. the temperature), and sample those extended configurations according 
to some stochastic rules. However, these methods asks for a prior distri- 
bution of the additional parameters (for example a temperature ladder in 
for parallel tempering method), which are usually estimated through some 
preliminary runs |19j . 

As noted by many authors, simulated annealing strategies can be used 
to compute exactly ratios of partition function (free energy differences), 
through an explicit computation of importance weights of nonequilibrium 
simulations, often referred to as Jarzynski's equality POII^ . 

We present here a complementary approach to the above simulated an- 
nealing type strategies. It is similar to the method of known as "pop- 
ulation Monte-Carlo", in which multiple replicas are used to represent the 
distribution under study. A weight is associated to each replica, and re- 
samplings are performed at discrete fixed times to avoid degeneracy of the 
weights. This methodology is widely used in the fields of Quantum Monte 
Carlo ^ 1^ or Bayesian Statistics, where it is referred to as Sequentiel 
Monte Carlo IHj. Note that in the probability and statistics fields, each 
simulation is called a 'walker' or 'particle'; we use here the name 'replica', 
which is more apppropriate to the Molecular Dynamics context. 

Our method extends the population Monte-Carlo method to the time- 
continuous case. It consists in running M replicas of the system in parallel, 
resorting typically to a stochastic dynamic, and considering exchanges be- 
tween them, according to a certain probability depending on the work done 
on each system. This procedure can be seen as automatic time continu- 
ous resampling, and all replicas have the same weight at any time of the 
simulation. This method drastically increases the number of significative 
transitions paths in nonequilibrium simulations. These heuristic explana- 
tions are precised in section |21 The set of all replicas (or walkers) is called 
an 'Interacting Particle System' (IPS) 10;, and can be seen as a genetic 
algorithm where the mutation step is the stochastic dynamics considered. 

The article is organized as follows. We first precise classical simulated 
annealing type methods in section We then describe the associated IPS 
method in section [2 as well as its numerical implementation. Possible ap- 
plications and some numerical results are then presented in section |31 and 0] 
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1 Simulated annealing type methods 



Consider a family of Hamiltonian functions Hx : T*A4 — > M indexed by a 
parameter A G [0, 1]. The corresponding Hamiltonian dynamics are 



dq 




'dt~ 


dp ' 


dp 


dHx 


~(E~ 


dq 



(5) 



The family (-ffA)AG[o,i] indexes a path between the original state described by 
a Hamiltonian Hq and the final state charaterized by a Hamiltonian Hi. A 
canonical probability measure ^\ can be associated to each Hamiltonian Hx : 

dfix{q,p) = ^e-''''^^'''PUqdp, (6) 
where the normalizing constant Zx is 

Zx= [ e-f^^^^^'P^qdp. 

Jm 

We will also denote by F{X) = —f3^^ InZ^ the corresponding free energy. 

We wish to sample according to dfii, from an initial, easily obtained 
sample of dfiQ. For each replica of the previous sample, the corresponding 
configuration of the system is brought slowly to the end state along a path 
(A(t))(g[o^T] foi' a time T > 0. The final sample of configurations is hopefully 
close to dni. 

Typically, we can consider Hx{q,p) = (1 — X)HQ{q,p) + XHi{q,p) (e.g. 
when performing a change temperature Ho{q,p) = H{q,p), Hi{q,p) = 
^ H{q,p)). It can also represent a modification of the potential, sometimes 
called 'alchemical transition' in the physics and chemistry litterature. The 
folding of a protein could be studied this way for example, by setting initially 
all the long-range interactions to zero, whereas the final state corresponds 
to a Hamiltonian were all interactions are set on. In this case, 

Hx{q,p) = ]^p^M-^p + Vx{q). 

The simulated annealing like strategies can also be extended to the re- 
action coordinate case j25j . In this case, the initial and the final state are 
indexed by some order parameter z{q). 

1.1 Markovian nonequilibrium simulations 

The usual way to achieve this method is to perform a time inhomogeneous 
irreducible Markovian dynamic 

t^xl^\ X^^-Mo, (7) 
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for t £ [0,T], and a smooth schedule t i— > A(t) verifying A(l) = and 
A(T) = 1, and such that for all given A G [0, 1], the Boltzmann distribution 
dfi\ is invariant under the dynamics t X^. 

The variable x can represent the whole degrees of freedom {q,p) of the 
system, or only the configuration part q. Depending on the context, the in- 
variant measure ^ will therefore be the canonical measure @, or its marginal 
with respect to the momenta, which reads 

d/iA((z) = ^e-^^-('?)(i<Z, (8) 

with 

Jm 

When we do not wish to precise further the dynamics, we simply call dfix{x) 
the invariant measure, and x the configuration of the system. The actual 
invariant measure should be clear from the context. 

For all t € [0,T], the dynamic (O will be usefully characterized by its 
infinitesimal generator Lx{t)^ defined on a domain of continuous bounded 
test functions by: 

L,(,)(v.)(x) = hm i(E(^(X,f = x) - ^{x)) 

The invariance of /^A(t) under the instantaneous dynamics can be expressed 
through the balance condition: 

V</5, /iA(t)(^A(t)(93)) = 0. (9) 
The dynamics we have in mind are (for a /ixed A G [0, 1]): 
• The hypo-elliptic Langevin dynamics on T*M 

dqi 
dpi 

where Wt denotes a standard 3A^-dimensional Brownian motion. The 
paradigm of Langevin dynamics is to introduce in the Newton equa- 
tions of motion © some fictitious brownian forces modelling fluc- 
tuations, balanced by viscous damping forces modelling dissipation. 
The parameters a, > represent the magnitude of the fluctuations 
and of the dissipation respectively, and are linked by the fluctuation- 
dissipation relation: 

a = {2i/pfl\ (11) 



dH 



dp 
dq 



iqt,pt)dt, 



{qi,pi) dt - iM-^p^ dt + a dWt, 



(10) 
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Therefore, there remains one adjustable parameter in the model. The 
infinitesimal generator is given by: 



• The elhptic overdamped Langevin dynamic^ in the configuration space Ml : 

dq^ = -VVxiq^) dt + a (12) 
where the magnitude of the random forcing is given here by 



" -Si- 

The corresponding infinitesimal generator is given by: 

Let us remark that the overdamped Langevin dynamic 1)12^ is obtained 
from the Langevin dynamic (|10j) by letting the mass matrix M go to 
zero and by setting ^ = 1, which amounts here to rescaling the time. 

It is well known that, for a fixed A G [0, 1], these dynamics are ergodic under 
mild assumptions on the potential V 

When the schedule is sufficiently slow, the dynamics is said quasi-static, 
and the law of the process X^^^^ is assumed to stay close to its local steady 
state throughout the transformation. As said before, this is out of reach at 
low temperature (more precisely, large deviation results p!B| ensure that the 
typical escape time from metastable states grows exponentially fast with P, 
which implies quasi-static transformations to being exponentially slow with 
P)- 

It is therefore interesting to consider approaches built on the simu- 
lated annealing formalism, but able to deal with reasonably fast transition 
schemes. 



1.2 Importance weights of nonequilibrium simulations. 

For a given nonequilibrium run X^^^^ we denote by 

Wt = l^^^{X^^'^)X'{s)ds (13) 

""^This dynamics is actually known as the 'Langevin dynamic' in the probability and 
statistics fields. We adopt here the physical names of these stochastic processes, which 
are more natural when dealing with molecular dynamics. 
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the out of equilibrium virtual work induced on the system on the time sched- 
ule [0, t] . The quantity Wt gives the importance weights of nonequilibrium 
simulations with respect to the target equilibrium distribution. Indeed, it 
was shown in j2Uj that 

]E(g-/3m) ^ g-/3(F(A(t))-F(0))_ ^^4^ 

This equality is known as the Jarzynski's equality, and can be derived 
through a Feynman-Kac formula JH]- Differentiating the un- normalized 
Boltzmann path 1 1— > = e~I^^Hi)^^^ dx with respect to t: 

and using the balance condition ©, 

dtU^^t^i^) = n;,(,) (^L^^t){^) - /3^^X'{t)^^ , (15) 



it follows 

,mi^)o.^^ = EUx^^%-^^A. (16) 

Therefore, taking ip = 1, 

]E(e-/3m) ^ g-/3(F(AW)-F(0))_ ^^7^ 

Jensen's inequality then gives 

E(m) > F{X{t)) - F(0). (18) 

This inequality is an equality if and only if the transformation is quasi-static 
on [0,t]; in this case the random variable Wt is actually constant and equal 
to AF. When the evolution is reversible, this means that equilibrium is 
maintained at all times. 

As an improvement, we will use the importance weights Wt to perform 
a selection between replicas. 



2 The Interacting Particle System method 

Our strategy is inspired by the re-sampling methods in Sequential Monte 
Carlo (SMC) literature [HI In this time continuous context, the idea 
is to replace importance weights of simulations performed in parallel, by a 
selection operation between replicas. 

We first describe the IPS approximation in section 12.11 as well as con- 
vergence results of the discretized measure to the target measure. A justifi- 
cation through a mean-field interpretation is then presented in section 
The numerical implementation of the IPS method is eventually discussed in 
section 12.31 
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2.1 The IPS ans its statistical properties 

Recall that the potential of mean force defined by 



is the average force applied to the system during an infinitely slow trans- 
formation. It can be used in a thermodynamic integration to compute free 
energy differences: 

F(1)-F(0)= r F^t)\'{t)dt. (19) 

The first step is to rewrite the Feynman-Kac formula (|16|) by introducing a 
dichotomy when a replica is receiving either excess or deficit work compared 
to the potential of mean force. 

To this end, we define respectively the excess and deficit force, and the 
excess and deficit work as 



11',-= / /f(X^W)A'(s)rf., = / /f (X.*M)A'W ds, (20) 

JQ Jo 

and rewrite 

We now present the particle interpretation of (|21() enabling a numerical 
computation through the use of empirical distributions. Consider M Marko- 
vian systems described by variables {0 < k < M). We approximate the 
virtual force by 

k=l 



and the Boltzmann distribution by 

M 

M 



1 



k=l 



which are their empirical versions. This naturally gives from definitions (|2()j) 
empirical approximations of excess/deficit forces j^^''^^/'^^ a,nd works 

The replicas evolve according to a branching process whith the following 
stochastic rules (see [HS] for further details) : 
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Process 1. Consider an initial distribution {Xq,...,Xq ) generated from 
(i/io(x). Generate idependent times t^'^jT^''^ from an exponential law of 
mean f3~^ (the upperscripts b and d refer to 'birth' and 'death' respectively) , 
and initialize the jump times T^^^ as 
For{)<t<T, 

• Between each jump time, evolve independently the replicas accord- 
ing to the dynamics ((TJ; 

• At random times T^'^^ defined by 

rpk,d rpk,d 'n+1' 

-'n + 1 -'" 

an index Z € {1, . . . , M} is picked at random, and the configuration of 
the k-th replica is replaced by the configuration of the l-th replica. A 
time ^•^ generated from an exponential law of mean 

• At random times T^'^^i defined by 

Ti^A;,de _ Ti/fc,de _ k,b 
rpk,d rj,k,d — 'n+1' 

-'n+l -'n 

an index I G {1, . . . , M} is picked at random, and the configuration of 
the l-th replica is replaced by the configuration of the k-th replica. A 
time generated from an exponential law of mean . 

The selection mechanism therefore favors rephcas which are samphng 
values of the virtual work Wt lower than the empirical average. The system 
of replicas is 'self-organizing' to keep closer to a quasi-static transformation. 

In \M\. several convergence results and statistical properties of the 
replicas distribution are proven. They are summarized in the following 
proposition: 

Proposition 2.1. Assume that {t,x) ^ ^^(x) is a continuous bounded 
function on [0, T] x T*Ai (or [0, T]x A4 in the case of overdamped Langevin 
dynamics), and that the dynamic ((T)) is Fellerian and irreducible. Then for 
all test function (p and any t G [0,T], 

• The estimator 

^*f,)Mexp (-/3^V*f,)A'(.)d.) (22) 

is an unbiased estimator of 

• the estimator iJ,^^^^{f) is an asymptotically normal estimator of iix(i^{f), 
with bias and variance of order M~^. 
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The proof follows from Lemma 3.20, Proposition 3.25 and Theorem 3.28 
of ^Uj (see also for further details). 

The unbiased estimation of un-normalised quantities is a very usual prop- 
erty in importance sampling type strategies. In discrete time SMC methods, 
it comes from the fact that at each time step, when operating re-sampling, 
each replica branches with a number of offsprings proportional in average to 
its importance weight. Unfortunately, there is no simple justification of this 
fact for time continuous IPS since the branching phenomenon is diluted in 
the continuity of time. 



2.2 Justification through a mean-field hmit 

In order to prove the consistency of the IPS approximation, we consider the 
ideal setting where the number of replicas goes to infinity (M +oo). This 
point of view is equivalent to a mean-field or Mc Kean interpretation of the 
IPS (denoted by the superscript 'mf'). In this limit, the behavior of any 
single replica, denoted by X™^, is then independent from any finite number 
of other ones. We shall consider the mean field distribution and force: 



Law(Xr') = 

^mf _ ,.mf f 9Hm \ 



The mean field excess/deficit force j™^'^^/^^ g^j^d works ^^^^'"^^^^^ are defined 
as in (|20|) . 

In view of Process ^ the stochastic process Xf^^ is a jump-diffusion 
process which evolves according to the following stochastic rules (some facts 
about pure Markov jump processes are recalled in the Appendix): 

Process 2. Generate from d^jLQ{x). Generate idependent clocks (r^, T^)n>i 
from an exponential law of mean j3~^ (the upperscripts b and d refer to 'birth' 
and 'death' respectively), and initialize the jump times T^l'^ as Tq = 0,Tq = 
0. 

ForO<t<T, 

• Between each jump time, t X™^ evolves according to the dynamic Q ; 

• At random times T^+i defined by 

the process jumps to a configuration x, chosen according to the proba- 
bility measure d^^^ {x); 
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• At random times T^+i defined by 

the process jumps to a configuration x, chosen according to the proba- 
bility measure proportional to f^^''^'^{x)dfi^(rpb ){x). 

Remark 1. Note that, in the treatment of the deficit work, we take in 
Process the point of view of the jumping replica; whereas in Process Q we 
take the point of view of the attracting replica which induces a branching. 

Prom the above probabilistic description, we can derive the Markov gen- 
erator of the mean-field process, given by the sum of a diffusion and a jump 
generator: 

where the jump generator ^mf is defined as 

J,^^^.{^){x)=/3X'{t) [ {^{y)-^{x)){ff^'''{x) + ff^''%y))^,f{dy). 
J M 

The non-linear evolution equation of the mean-field distribution /i™^ is then 

which gives, by a straightforward integration, 

dtl^fif) = l^f + /? (^Tf - A'(t)/ 

so that finally 

Since the un-normalized Boltzmann distribution is the unique solution of 
the above linear equation (seell5|). we obtain the following identities: 

,,mf _ ,, -r-mf _ -r- „mf,cx/de _ „cx/de /„q\ 

IJ-t — l^\{t)^ -^t — •^A{t); It ~ hit) ■ ^'^'^> 

This proves the consistency of the IPS approximation scheme. 
2.3 Numerical implementation 

In the previous section, we discretized the measure by considering an em- 
pirical approximation. For a numerical implementation to be tractable, it 
remains to discretize the time evolution. Notice already that the IPS method 
induces no extra computation of the forces, and is therefore unexpensive to 
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implement. However, although the IPS can be parallelized, the processors 
have to exchange informations at the end of each time step, which can slow 
down the simulation. 

There are several ways to discretize the dynamics ()lUp or H12|) . Some 
common schemes used in molecular dynamics are the Euler-Maruyama dis- 
cretization for (|12|) . and the BBK scheme |3] for We refer to 31 
alternative approaches in the field of molecular dynamics. In the sequel, we 
will denote by x^''' a numerical approximation of a realization of Xf^^, with 
X = q OY X = {q,p) depending on the context. 



Euler discretization of the overdamped Langevin dynamics. The 

Euler-Maruyama numerical scheme associated to (|12|) reads, when taking 
integration time steps At, 



qu+i =qn_^^ VV^iq^) + ^ — R^, (24) 

where (ii")„gN is a sequence of independent and identically distributed 
(i.i.d.) 3A^-dimensional standard Gaussian random vectors. The numeri- 
cal convergence of this scheme can be ensured in some circonstances 



Discretization of the Langevin dynamics. When considering an in- 
tegration time step At, the BBK discretization of (jlUj) reads component- 
wise (recall that the underscripts j refer here to the components of a given 
x^ = x= {q,p)), 



n+l/2 
Pj = 


n 


(-Vg,F((? 


< qr = 


q^ + At 


n+1/2 
Pj 

rrij 


< 


1 


I n+1/2 
[Pj - 


^ 2m, 



.p1 



At ^ 



(25) 



At 



2 V,^,y((?"+i)+a 




where the random forcing terms i?" (j G {1,...,A^} is the label of the 
particles, n is the iteration index) are standard i.i.d. Gaussian random 
variables. The fluctuation/dissipation relation (jllj) must be corrected so 
that the kinetic temperature is correct in the simulations 1^. To this end, 
we set 

Notice that the relation (jllj) is recovered in the limit At 0. 
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Discretization of the selection operation We consider for example the 
foUowing discretization of the force exerted on the k-th repUca on the time 
interval [iAt, (i + l)At]: 



dX 2\ dX ' ' dX 
The mean force is then approximated by 

j-M,At _ _j_ -^^ + l/2 

k=l 

To get a time dicretization of the IPS, Process ^ is mimicked using the 
following rules: 

• the time integrals are changed into sums; 

• the selection times are defined as the first discrete times exceeding the 
exponential clocks t^/'^. 

Further details about the numerical implementation can be found in j32j . 
Note that one can find more elaborate methods of discretization of the IPS 
(see j33|). but this one seems to be sufficient in view of the intrinsic errors 
introduced by the discretisation of the dynamics. 



3 Applications of the IPS method 

3.1 Computation of canonical averages 

The most obvious application of the IPS method is the computation of 
phase-space integrals, since an unweighted sample of the target Boltzmann 
distribution /xi is generated. The sample obtained can of course be improved 
by some additional sampling process (according to a dynamics leaving the 
target canonical measure invariant). This will decorrelate the replicas and 
may increase the quality of the sample. We refer to ^ for further precisions 
on sampling methods, to for some numerical tests assessing the quality 
of the sample generated as a function the transition time T and the number 
of replicas M, and to j32j for an application to pentane. 

3.2 Estimation of free energy difference. 

The free energy of a system is defined as 

F = -p-HnZ 

where Z is the partition function Z = j exp{—PH{q,p)) dqdp. It cannot be 
computed with a single sample of nx. Only free energy differences can be 
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computed easily. Since the free energy of certain states is known (This is the 
case for perfect gases, or for sohds at low temperature 1301)) the free energy 
of any state can in principle be obtained by an integration between a state 
whose free energy is known, and the state of interest. Usual methods to 
this end are Thermodynamic integration j2l], the free energy perturbation 
method and the related Umbrella sampling technique |3Z1) or Jarzynski's 
non equilibrium dynamics (also called 'fast growth') PUJ. It is still a matter 
of debate which method is the most efficient. While some results show that 
fast growth methods can be competitive in some situations jTB], other studies 
disagree |29| . However, a fair comparison is difficult since the dynamics used 
may differ (in Hamiltonian dynamics are used during the switching 
process), and more efficient fast growth methods techniques (using, e.g. path 
sampling 36^ JSI) are still under investigation. 

In the work of Jarzynski [20], M independent realizations {X^ , X^^) 
of a bare out of equilibrium dynamic ((Tj) are used to compute free energy 
differences through (|14j) . with the estimator 



An alternative (yet similar) estimator relying on the thermodynamical inte- 
gration^! which is also considered in jHS] (and is implicit in |2Uj ) is 



However, both estimators AFj and AFj suffer from the fact that only a 
few values of are really important. Indeed, because of the exponential 
weighting, only the lower tail of the work distribution is taken into account. 
The quality of the estimation then relies on those rare values, which may be 
a problem in practice (see e.g. EHI). 

In the case of interacting replicas, we use similarly 





where 





which shares by Proposition 12.11 the same statistical properties as AFj: 
AFjPS is asymptotically normal with bias and variance of order M~^, and 
the estimator q-P^^ips \g unbiased estimator of e"'^'^'^. 
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Our numerical comparisons often turned out to give similar free energy 
estimations for the IPS method and the standard Jarzynski method. How- 
ever, we have mostly considered the issue of pure energetic barriers, where 
the difficulty of sampling comes from overcoming a single high barrier. The 
observed numerical equivalence may be explained by the fact that the se- 
lection mechanism in the IPS method does not really help to explore those 
regions of high potential energy. 

However, when the sampling difficulties also come from barriers of more 
entropic nature {e.g. a succession of very many transition states separated 
by low energy barriers), the IPS may improve the estimation. Indeed, the 
selection mechanism helps keeping a statistical amount of replica in the 
areas of high probability with respect to the local Boltzmann distribution 
throughout the switching process (see the numerical example in section HTTj) . 
This relaxation property may be crucial to ensure at each time a meaningful 
exploration ability. These points are still under investigation. 

3.3 Initial guesses for path sampling 

The problem of free energy estimation is deeply linked with the problem of 
sampling meaningful transition path. In the IPS method, one can associate 
to each replica a genealogical continuous path (^s '^^'^)sg[o,t]- The latter 
is constructed recursively as follows for a replica k (for < t < T): 

• at each time t, set X^'^'^'^ = Xj^; 

• at each random time r„ when the replica jumps and adopts a new 
configuration (say of replica set (X^ '^°'')[o,t„] = (-'^l' ^™)[o,T„] • 

This path represents the ancestor line of the replica, and is composed of 
the past paths selected for their low work values. The study of the set of 
genealogical paths lies outside the scope of this article (see jH] for a discussion 
in the discrete time case). However, let us mention that for a given t £ [0, T], 
the set of genealogical paths is sampled, in the limit M — > oo, according to 
the law of the non-equilibrium paths {Xs^^^)s^[o,t] weighted by the factor 
^-l3Wt (^^ith statistical properties analogous to those of proposition 12. 
These paths are thus typical among non-equilibrium dynamics of those with 
non-degenerate work. Therefore, they might be fruitfully used as non-trivial 
initial conditions for more specialized path sampling techniques (as e.g. |38j ) . 

4 Numerical simulations 

4.1 A toy example 

Consider the following family of Hamiltonians (-ffA)Ae[o,i]- 

Hx{x) = y + AQi(x) + yQ2(x) + ^Q'six) + ^Q^ix) (27) 
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4/5 



Figure 1: Plot of some Hamiltonian functions, as defined by H27|). 



with 



-1 



8x2 + 1' 



Q2{X) 



-18 



8(x- 1)2 + 1' 
-84 



32(2; -3/2)2 + 1' ' 64(x- 7/4)2 + 1' 

Some of those functions are plotted in Figure ^ This toy one-dimensional 
model is reminiscent of the typical difficulties encountered when /^q is very 
different from /xi. Notice indeed that several transitional metastable states 
(denoted A and B in Figure ^) occur in the canonical distribution when 
going from A = to A = 1. The probability of presence in the basins of 
attraction of the main stable states of Hi (C and D in Figure ^ is only 
effective when A is close to 1. 

Simulations were performed at /? = 13 with the overdamped Langevin 
dynamics (|12() . and the above Hamiltonian family ()27() . The number of 
replicas was M = 1000, the time step At = 0.003, and A is considered to 
be linear: A(t) = t/T. Figure HI presents the distribution of replicas during 
a slow out of equilibrium plain dynamic: T = 30. Figure El presents the 
distribution of replicas during a faster dynamics with interaction: T = 15. 

When performing a plain out of equilibrium dynamics (even 'slow') from 
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D 




>L=4/5 ?L=1 



Figure 2: Empirical densities (in dots) obtained using independant replicas. 




Figure 3: Empirical densities (in dots) obtained using interacting replicas. 

A = to A = 1, almost all replicas are trapped by the energy barrier 
of these transitional metastable states (see Figure [2). In the end, a very 
small (almost null) proportion of replicas have performed interesting paths 
associated with low values of virtual work W. When using H16() to compute 
thermodynamical quantities, these replicas bear almost all the weight of the 
degenerate sample, in view of the exponential weighting. The quality of the 
result therefore depends crucially on these rare values. 

On the contrary, in the interacting version, the replicas can perform 
jumps in the configuration space thanks to the selection mechanism, and go 
from one metastable basin to another. In our example, as new transition 
states appear, only few clever replicas are necessary to attract the others in 
good areas (see Figure|31). In the end, all replicas have the same weight, and 
the sample is not degenerate. Notice also that the final empirical distribution 
is fairly close to the theoretical one. 

We have also made a numerical estimation of the error of the free energy 
estimation, with 40 realizations of the above simulation. The results are 
presented in Tabled and show an important reduction of standard deviation 
and bias up to a factor 2 when using the IPS method. 
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Method 


Bias 


Variance 


Plain 


+0.25 


0.19 


Interacting 


+0.15 


0.10 



Table 1: Error in free energy estimation. 



4.2 Gradual Widom insertion 

We present here an application to the computation of the chemical potential 
of a soft sphere fluid. This example was considered in jl6l I29j for example. 
We consider a two-dimensional (2D) fluid of volume \^l\, simulated with 
periodic boundary conditions, and formed of particles interacting via 
a pairwise potential V. The chemical potential is defined, in the NVT 
ensemble, as 

dF , , 

where F is the free-energy of the system. Actually, the kinetic part of the 
partition function Z can be straightforwardly computed, and accounts for 
the ideal gas contribution In the large limit, the chemical potential 
can be rewritten as ^3] 

whith 

/^id = In 



(Ar + 1)A3 

where A is the "thermal de Broglie wavelength" A = /i(27rm/3~^)~^/^ (with 
h Planck's constant). The excess part //ex is 

« - / Wexp(-/3y(g^^^))d,^^n 

^ '\ |17|^exp(-/3y(g^))dg^ J' ^'^^ 

where 1^((?^) is the potential energy of a fluid composed of A^ particles. We 
restrict ourselves to pairwise interactions, with an interaction potential ^. 
Then, V((?^) = Y.l<^<■J<N^{W^-(^3\)■ Setting 7r(g^) = exp(-/3y(g^)) and 
Ay(g^,g) = 1/(g^+i) - V{q^) with g^+i = (g^,g), it follows 



/.ex = -a-' In ( ^ / e-^^^(«'^")d7r(g^) dq ) . (30) 



f7 



The formula (jHOJ can be used to compute the value of chemical poten- 
tial using stochastic methods such as the free energy perturbation (FEP) 
method j39|. In this case, we first generate a sample of configurations of the 
system according to vr, and then evaluate the integration in the remaining 
q variable by drawing positions q of the remaining variable uniformly in fi. 
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Another possibility is to use fast growth methods, resorting to the fol- 
lowing parametrization 

N+l 2 ^+1 2 

1=1 i=l 

(31) 

In this case, the interactions of the remaining particle are progressively 
tmrned on. 

As in |16[ I29j , we use a smoothed Lennard- Jones potential in order to 
avoid the singularity at the origin (Let us however note that, once the par- 
ticle is inserted, it is still possible to change all the potentials to Lennard- 
Jones potentials, and compute the correponding free-energy difference). The 
Lennard Jones potential reads 

and the modified potential is 

{a-6r2, 0<r<0.8cj, 
<^LJ (r) + c{r - rc) - d, 0.8^7 < r < r^, 
0, r > Tc. 

The values a, b, c are chosen so that the potential is C^. The distance Vc is 
a prescribed cut-off radius. We considered the insertion of a particle in a 
2D fluid of 25 particles, at a density pa^ = 0.8, with Vc = 2.5 a, /3e = 1, 
At = 0.0005, and a schedule A(t) = t/T where T is the transition time. The 
results are presented in Table |5J for different transitions times, but at a fixed 
computational cost, since MT is constant. Some work distributions are also 
depicted in Figure 01 A reference value was computed using FEP, with 10® 
insertions, done by running M = independent Langevins dynamics for 
the system composed of N particles, for a time tpEP = 50 (after an initial 
thermalization time to decorrelate the systems) , and inserting one particle at 
random after each time-step. The reference value obtained is = 1-32 k^T 
(±0.01 kBT). 

As can be seen from the results in Table 121 the IPS algorithm has a com- 
parable accuracy to Jarzynski's estimates. However, the work distribution 
is very different, and has a gaussian shape for all switching rates considered, 
whereas the work distribution obtained through the fast growth method are 
much wider (see in particular Figure |31 (Left)), so that the relevant part of 
the work distribution (the lower tail) is only of small relative importance. 
Of course, the width of this distribution decreases as the transition is made 
slower. 
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Method 


M = 10^ 
T = 1 


M = 5 X 10"^ 
T = 2 


M = 2 X 10"^ 
T = 5 


M = 10* 
T = 10 


SA 


1.32 


1.29 


1.33 


1.36 


IPS 


1.37 


1.35 


1.29 


1.34 



Table 2: Free energy estimation for one realisation of each method, depend- 
ing on the switching time T and the number of replicas M used, keeping 
MT constant. The reference value obtained through FEP is /Xex = 1-31 k-^T 
(±0.01 k^T). Notice that the results are quite comparable. 




Figure 4: Left: Comparison of the work distribution for T = 1. Right: 
Comparison of the work distributions for T = 10. The IPS results appear 
in darker colors. Notice that the bare simulated annealing process work 
distributions come closer to a gaussian shape, and has a reduced variance 
as T increases, as expected. 
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Appendix : Pure jump processes 

Consider a Markov process Xt of infinitesimal generator 



where for each time t, at is a bounded positive function and iit a probability 
measure. Denote by {Tn)n>i the jump times (with Tq = 0), and (r„)„>i 
independant clocks of exponential law of mean 1. Then the system evolves 
according to the following stochastic rules: 

• the jump times are defined by 



• at jump times, the process jumps to a configuration chosen according 
to the probability measure (i/^rn+i(y)- 





n 
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